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Shape Memory Alloys (SMA) are materials with increasing range of applications in engineering, 
aerosapce, and biomedical industries. They possess unique properties of being able to recover their 
original shape after permanent deformations. These materials can directly transduce thermal en- 
ergy into mechanical and vice versa. The key to "pseudo-elastic" behaviour of these materials and 
the "shape memory effect" [8,26] is held by the first-order martensitic phase transformations in 
these materials. Indeed, drastic changes in their properties are originated from their microstruc- 
ture or phase combinations. Via the phase transformation, the microstructure of the material can 
be switched among various combinations between austenite, martensite variants, or their mixtures. 
Austenitic phase is a more symmetric phase of the crystallic lattice, prevailing at high temperature, 
while martensite is a less symmetric, low-temperature phase [1,8,26]. Between these two critical 
situations, austenite and martensite might co-exist providing a typical example of phase combi- 
nations. Upon external loading, one phase combination can be switched to another. If the mate- 
rial is constrained at the boundary, a specific ("self-accommodating") combination of different 
phases will be established such that the bulk energy in the constrained domain will be minimized 



The mathematical framework for modeling phase combinations in shape memory materials is 
based on the solution of the variational problem with respect to a frame-indifferent non-convex 
free energy function 9) ([5,6,11,19,18] and references therein) : 



where Q is the reference configuration associated with the considered material, Y is the deforma- 
tion tensor, and 9 is the temperature of the material. Hence, by minimizing from (1), we 
minimize the bulk energy of the considered structure, as a functional of the deformation tensor, at 
temperature 9. This procedure, performed particularly often at the mesoscale [7,10,11,19,18], has 
several known difficulties. It is known that in the general case the variational problem given by 
Eq.(l) may have infinitely many minimizers ([7,19,10] and references therein). On the other hand, 
a more precise definition of the free energy that allow us to account for interfacial energy effects 
by introducing gradients of the order parameters is a highly non-trivial task ([1,18] and references 
therein) connected with additional difficulties. The problem can also be regularized by assigning 
simplified, e.g. affine, boundary conditions. For example, for the simulation of simple laminated 
microstructure, the affine boundary conditions are constructed by assuming that the deformation 
gradient of the material on the domain boundary is a linear combination of its equilibrium defor- 
mation gradients: 



[1,18,19]. 




(1) 



Y{x) = (AFo + (1 - X)RFi)x, x e dQ, 



(2) 
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where i? is a rotation matrix satisfying a twinning equation (see [7, 1 0, 1 1 , 1 9]) and A is the thickness 

of layers in the laminated microstructure. Fq and Fi are the deformation gradients that minimize 
the local energy function (for square to rectangular transformation) when the deformation gradient 
VF takes any of the following following values: 




r- 

where is a local minimum of the free energy function (f){Y, 9). In the general case, however, the 
above affine boundary conditions may not be appropriate. 

In solving problem (1), we have to face also numerical challenges connected with several local 
minima, resulted from phase mixtures under low and moderate temperature regimes, and non- 
convexity of the problem ([4,10,19,30] and references therein). Minimization procedures based 
on conventional local search methods with randomly chosen initial guesses [9] may not lead to 
a satisfactory result in those cases where the solution space is discretized with a large number of 
node points. 

In the present paper, we construct a mathematical model for the simulation of phase combinations 
in SMA materials on the basis of the Landau theory. We associate the phase combination with 
the global minimizer of the bulk energy of the SMA structure with the prescribed mechanical 
boundary conditions. We develop a hybrid optimization strategy consisting of two main steps. 
Firstly, we apply the Genetic Algorithm (GA) and its global exploration capacity to obtain an 
initial estimation of the global minimizer. Then, we apply the quasi-Newton method to refine 
such an estimation locally. Finally, the developed procedure is demonstrated by several numerical 
examples simulating phase combinations in SMA materials. 



2 Landau Free Energy Function and Variational Formulation 



For the modeling of phase combinations in SMAs structures, the first task is to characterize dif- 
ferent phases, which is different from one material to another. Here our mathematical model will 
be based on the square to rectangular transformations. In Fig.l (a) we give a schematic represen- 
tation of this case where the square lattice is the austenite, while two rectangles are the martensite 
variants. The square to rectangular transformation could be regarded as a 2D analog of the cu- 
bic to tetragonal or tetragonal to orthorhombic transformations observed in general 3D cases in 
NbsSn, InTl, FePd alloys and some copper based SMAs ([17,18] and reference therein). The 
analysis of this 2D transformation is a first step in understanding more complex cubic to tetrago- 
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nal and tetragonal to orthorhombic transformations. Most numerical studies of the dynamics of the 

phase transitions up to date have been concentrated on the analysis of the formation and growth 
of the micro structure ([1,17,18] and references therein). Such studies have been focused on the 
mesoscale under either periodic or affine boundary conditions. 

In what follows, we base our consideration on the Landau theory of phase transformation. Ac- 
cording to this theory, the basis of any nonlinear continuum thermodynamical model for phase 
transformation is a non-convex /ree energy function [23,17,18]. The local minima of the free en- 
ergy function with respect to the strain tensor (or deformation gradients) correspond to the stable 
and mesostable state at a given temperature, while the microstructure in the domain of interest can 
be described by the minimizer of the bulk energy in the domain. One of the simplest realization 
of this idea in the context of SMAs is based on the Helmholtz free energy ^ ([13,23,18,22]and 
references therein): 



where iJjo{9) models thermal field contributions, tpi{9)ip2{£) models shape memory contributions 
and ilJ3{e) models mechanical field contributions, e = du/dx (u is the displacement) is the strain 
which is chosen as the only order parameter in the ID case. The thermal field contributions ipo can 
often be modelled as follows [13,28] : 



where is the specific heat constant. 

The numerical analysis of a system of conservation laws based on the above representation of 
the free energy function (and a more general one) has been recently reported in detail in [22]. A 
conservative numerical scheme was constructed for the solution of the problem. It was noted that a 
standard energy inequality technique, applied to the convergence analysis of the scheme, can lead 
to quite restrictive assumptions. In [22] it was shown how such assumptions can be removed. 

This work focuses on the practical development of an algorithm suitable for the simulation of 
phase combinations. In this context we note that a free elastic energy functional (denoted further 
by F), similar to the one discussed above, was established earlier to characterize the austenite at 
high temperature and the martensite variants at low temperature in SMA patches, specifically for 
the square to rectangular transformation where the Landau free energy function Fi was modified 
[17,18,1,28]. Recall that for the square to rectangular transformation we have to deal only with two 
martensite variants and only one order parameter [13,17,18] in order to characterize the martensite 
variants and austenite in a 2D domain. Following previous works on the subject ([13,17,18,28] and 



(4) 



■00 = —CyOlnO, 



(5) 
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(6) 



^ i=i ^ i=i 



where V is the gradient operator, ^2,0^ i = I, 3,2, and 0^3 are material-specific coefficients, 

and ei, 62, 63 are dilatational, deviatoric, and shear components of the strains, respectively, defined 
as follows: 

ei = (7711 + 7722) /v^, 



where Ui is the displacement in the i*" direction in the Cartesian system of coordinates, jc= 
(xi, X2, X3) are the coordinates of a material point in the domain of interest. It is known that in 
this formulation the deviatoric strain 62 can be chosen as the order parameter. More precisely, 62 
and 63 are two-component order parameter strains that have been discussed before in [1,2,3]. 

In the above formulation the Ginzburg term Fg is the term proportional to square of strain gradi- 
ents. This term is often included to account for the presence of domain walls. It is essential in sim- 
ulating phase growth and several other phenomena ([17,18]and references there in), but it can be 
ignored if we are interested only in the macroscopic phase combinations of the SMAs patch under 
mechanical loadings (indeed, the simulation scale in this case is too coarse to capture mesoscale 
structures). Furthermore, the energy contribution of this term is typically small compared to other 
terms. 




(7) 



(^12 + ^21) /2. 



The Cauchy-Lagrangian strain tensor 77 is given by its components 




(8) 
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In order to be able to model the entire range of different phase combinations that correspond to 
different temperatures, the material parameter A2 is assumed to be temperature dependent A2 = 
0-2(0 — 00)7 where Oq is a critical temperature, responsible for the appearance of an addtional 
minimum and corresponding to the austenitic phase when temperature increases. In Fig.l (b) we 
present the plots of the Landau free energy function, defined in this way, for the entire range of 
temperatures of interest (the material Au23Cu3oZn47). We observe that the function has two local 
minima at low temperatures (210°), which correspond to two (rectangular in the interpretation of 
Fig.l (a)) martensite variants, while only one minimum at the center corresponds to the (square) 
austenite phase when the temperature is high (270"). When the temperature is in between two 
critical values (e.g., at around 245" in the figure), we observe that there are three local minima 
demonstrating co-existence of metastable and stable phases. 

Now, if we take the thermal contribution ^0 the same as in the ID case, the final form of the 
Helmholtz free energy function for the square to rectangular transformation will take the following 
form: 

^{9,8) = -c^eine + ^el + + F„ F, = ^ {9 - 9o) el - + |e^. (9) 

By substituting the free energy function into the variational problem given by Eq.(l), the phase 
combination problem can be written as the following variational problem. Given temperature 9, 
find the displacements Ux and Uy (in the x and y direction) that minimize the bulk energy: 

W{ux, Uy) = J (-c„e In 9 + ^el + ^ej + ^{9- 9^) el - ^et + |e^) dV. (10) 

Under a given temperature, the contribution of the thermal field will not change the profile of 
the local free energy function, but rather only shift it upwards or downwards. Hence, if the applied 
external force in 2D is given by its components f^, fy, the final problem to solve can be formulated 
as follows: 

J \ 2^^ y ^ y ^ "'^ ^2 - + — - fxUx - fyUy j dV ^ mm . 
n 

The above model is reduced to the well-known Falk model in the ID case [13] that can be written 
with respect the only order parameter e = du/dx: 

W{u) = / (I - ^0) - ^e' + ^e' - /«) dV, (12) 
n 
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As discussed in Section 1, we supplement the model (1 1) by appropriate boundary conditions. For 
all examples discussed in Section 4, these are clamped boundary conditions: 

u^ = Uy = 0, atx = xi, Xr, ory = yt, yb, (13) 



where xi and Xr are the left and right boundaries along the x direction, yt and y^ are the top 
and bottom boundaries along the y direction, as sketched in Fig. 3. External forces vary and are 
specified in Section 4. 



3 Numerical Implementation Based on Hybrid Optimization 



The above variational problem is non-convex and its solution can only be obtained by numerical 
methods. The procedure developed in this section consists of two main steps: firstly, the variational 
problem is converted into a nonlinear minimization problem by spatial discretization, and then the 
solution of the resulting problem is sought by a hybrid optimization strategy. 



3. 1 Spatial Discretization Procedure 



For the spatial discretization, we employ the Chebyshev pseudospectral approximation (e.g., [4]) 
on a set of 2D Chebyshev points (xj, yj) in the 2D domain of interest Vt = [—1, 1] x [—1, 1]: 

Xi = cos{—),yj = cos{—), i,j = 0,1,... ,N, (14) 



where N + 1 is the total number of nodes in one direction. Other structures of interest can be 
mapped onto the domain Q by a linear transformation. Using the constructed grid, the displace- 
ments in the patch can be approximated as follows: 

N N 
i=0 j=0 



where f{x, y) represents either function or Uy, fij is the function value at (xj, y^). Functions 
^i{x) and ^j{y) are the i*^ and Lagrange interpolating polynomials along the x and y directions 
respectively. 
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Having obtained (15), the derivatives of function /, df{x, y) / dx and df{x, y) / dy, can be obtained 
by calculating d^i{x)/dx and d^j{y)/dy. Following the standard technique found, e.g., in [27,4], 
all the differentiation operators in Eq.(l 1) (or Eq. 12) can be written in the matrix forms: 

F, = D,F, Fy = DyF, (16) 



where F^ and F are vectors collecting all values of the derivative df/dx and the function / 

at {xi,yj), respectively, and similarly for Fy. The differentiation matrices Dj. and Dy can be 
calculated using the approximation given by Eq.(15). For instance, the differentiation matrix 
for the Falk model takes the following form: 



r 2Ar2 + 1 

6 

2N'^ + 1 
6 

"2(1 -,x|) 

. Cj {Xi — Xj) 



J 

3 



0, 

1,2, 



(17) 



N-l, 



i^j, i,j = l,2,...,A^-l, 



where Ci — 2 for: i — Q,N and Cj = 1 otherwise. Such matrices have dimensionality (iV + 1) x 
(AT + l). 

The bulk energy is given by an integral operator with the local free energy function as its inte- 
grand. We use the same set of points as chosen for the derivative approximation for constructing a 
numerical integration formula for integral. In particular, we use Chebyshev collocation nodes and 
the resulting quadrature formula is constructed by using the Chebyshev-Lobatto rule [4,27]. For 
example, the formula for integration in the x-direction, generically represented by 

/ f{x)dx ^J2^ifi^i)^ (18) 



is exact for any polynomials with an order less than 2N — 1. In (18) weight coefficients Wi are 
defined in the standard manner [4,15]. 

By substituting the approximations of all the differentiation operators and integral operators into 
the variational problem, we convert the original problem into the following minimization problem: 



find ulj, to minimize: W{ulj, ul^), 



(19) 
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where ujj and uf j stand for the values of Ux and Uy at node (xi^yj), respectively. Note that 

W{u} ,j, v!fj) is just the discretized bulk energy and is a nonlinear algebraic function of u} .j and 
ufj. The resulting problem has 2x(A^ + l)x(A^+l) variables in total. Given the prescribed 
boundary conditions, we have 2 x (iV — 1) x (iV — 1) variables in total. 

The problem to be solved is non-convex minimization problem with strong nonlinearity. It is not 

easily amenable to conventional gradient-based minimization methodologies due to multiple local 
minima. On the other hand, the genetic algorithms (GA) can be helpful in locating an approxi- 
mation to the global minimum, giving an initial approximation to gradient-based procedures. In 
what follows, we combine these two ideas by employing a hybrid optimization strategy that takes 
advantage of the global exploring capability of the GA and the local refinement accuracy of the 
quasi-Newton method for non-convex problems. 



3. 2 Genetic Algorithm Locates Initial Approximations 

The GA is a well established methodology for global optimization ([16,20,14] and references 
therein). We highlight here only its main features in the context of our problem. The GA maintains 
a population of individuals (chromosomes), say P{n), for generation n and each chromosome con- 
sists of a set of genes, where each gene stands for a parameter to be estimated. One chromosome 
represents one potential solution to the minimization problem. Each chromosome is evaluated to 
give some measure of its fitness according to the bulk energy defined in the previous section. 
Some chromosomes undergo stochastic transformations by means of genetic operations to form 
new chromosomes. Recall that there are two transformations in the GA: crossover, which creates 
new chromosome by combining parts from two chromosomes, and mutation, which creates a new 
chromosome by making changes in a single chromosome. New chromosomes, called offsprings 
S{ri), are then evaluated. A new population is formed by selecting fitter chromosomes from the 
parent population and the offspring population. After some generations, the algorithm converges to 
the fittest chromosome, which represents an estimated optimal solution to the problem [16,20,14]. 

Generation of initial chromosomes: In most of the GAs, the chromosomes are created by randomly 
choosing the genes in a given range. For the current problem, this is not an effective way. Indeed, 
if the strain in solid structures under consideration is assumed to be not very large, while displace- 
ments may vary slowly and be represented by smooth functions, there is no good reason to include 
high frequency oscillations in displacements. At the same time, if the displacement values are ran- 
domly chosen on our discrete grid, high frequency oscillations may well be pronounced. Hence, 
we smooth the randomly chosen chromosomes by the following filter operation: 



Us — IdlsUri 



(20) 
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where Ur is the vector collecting randomly chosen displacement values at all discretization nodes, 
while Us is the resultant smoother profile of the displacements after the filter operation. The 
matrix Ig is an interpolation matrix which maps the given data at the discretization nodes onto 
a sparser grid, using the least square approximation, while the matrix 1^ is another interpolation 
matrix which maps the data on the sparse grid back to the given discretization nodes. These two 
matrices can be constructed by using the Chebyshev collocation method, as discussed above. The 
product of these two interpolation matrices could be regarded as a filter to remove high frequency 
components. 

Crossover. Crossover operator in the GA is for producing new children chromosomes from chosen 
parent chromosomes. For a given pair of parent chromosomes Xi and X2, the offspring is obtained 
by the so called Intermediate recombination [20,15]: 

Oi = Xia + X2{1 - cx) , (21) 



where a is a vector with the same size as x, and all its entries are randomly chosen independently. 
This operation is capable of producing variables slightly larger than the hypercube in the solution 
space defined by the parents, but confined by the parameter ex. A typical range for the parameters 
in a is [—0.25, 1.25], which is used in the current paper. 

Mutation: The mutation operation is randomly applied with a low probability, typically in the 
range 0.001 to 0.01, and modifies genes in the chosen chromosome. Its role in the GA is to make 
sure that the probability of searching any potential solution (any points in the solution space) is 
nonzero, and is often regarded as a measure to recover good genetic material that might be lost 
through the operation of selection and crossover [20,16,14,15]. 

In practice, the mutation operation is carried out by replacing a randomly chosen gene by a ran- 
domly chosen new value in the given range. The mutated chromosome itself is randomly chosen 
from the current population with a given probability. In the current problem, the mutation op- 
eration is applied to one chromosome in one generation, which means the probability for each 
chromosome is one divided by the number of chromosomes. 

Generation-alteration method: Generation alteration method determines how to evolve the current 
generation to the next, which means how to select pairs of parents for producing children by 
crossover and mutation operator, and how to select parents in the currents population that survive 
in the next generation. 

Taking into account their fitness, it is obvious that each chromosome should have a probability 
determined by the associated function value. Here a simple linear map method is used as follows: 
rank all the chromosomes in the current generation in order of decreasing function values, then the 
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probability of being chosen for a specific chromosome i is calculated as: 



I 



(22) 



Pi = 



2{N + 1)N' 



where i is the position of the chromosome in the rank, while N is the number of chromosomes 
in the rank. Then the parents for survival and crossover are chosen using the above calculated 
probability. While all the chromosomes have nonzero probability to be chosen, fitter chromosomes 
have a better chance. 

The results of application of this procedure to the simulation of SMA phase combinations are 
reported in the next section. 



3.3 Quasi-Newton Method Refines the Solution 

The output of the GA is used as the initial guess for local search methods. In what follows we apply 
the quasi-Newton method to the minimization of the bulk energy given in Eq. 1 1 . Let's denote any 
potential minimizer by x. The bulk energy W{x) is minimized by x* if it satisfies VW{x*) = 0. 
To achieve this, we can organize the following iterative process. Let at the general /c*'^ step of the 
iteration, the potential solution is x^. Then, the task is to estimate the next Xk+i = Xk + dk such 
that VW{xk+i) = which means that 



where q{xk) = VW{xk) is the gradient of the bulk energy at Xk, Vq{xk) = V'^W{xk) is the 
Hessian matrix, while dk is the search direction. To make VW{xk+i) — 0, the search direction 
should satisfy: 



This standard Newton-type procedure is computationally expensive and has other well-known 
drawbacks [4,30] that preclude us from using it in the context of our problem. Instead, we apply 
the quasi-Newton procedure by constructing a matrix (at the A;*^ iteration), an approximation 
to the Hessian matrix that satisfies the following condition: 



VW{xk + dfe) = q{xk + </fe) = q{xk) + Vq{xk)dk = 0, 



(23) 



dk = -(Vq(xfc)) q{xk). 



(24) 



q{xk) + Bkdk = 0. 



(25) 
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Provided with the initial guess from the GA, the local search method using the quasi-Newton 
method is organized in a standard manner: at the general k^^ iteration (Xk is the current estimated 
solution): 

(a) Find a descent direction d!^ using Bk by the following formula: 

= -B^^VW{x''), (26) 

(b) Compute the acceleration parameter ak by line search, find ak such that W{x'' + akd!') is 
minimized. 

(c) Update the potential solution: 

Finally, the update of the approximation of the Hessian matrix B^+i is organized by using the 
BFGS update (Broyden-Fletcher-Goldfarb-Shanno update, e.g., [4,30,21]): 

tlk+1 - -Dfe ,tr . ^ 'TTrT' ^^^-^ 



where x/i^ = VW{xk+i) — VW (x^) . To initiate the iteration process, Bq for the first step is taken as 
the identity matrix that corresponds to the deepest descent methodology. For the current problem, 
we did not apply the line search for the acceleration parameter a^. Instead, a small value, regarded 
as a relaxation factor, was assigned to a^. 



4 Numerical Examples: Phase Combinations in SMA with Hybrid Optimization Procedure 



By combining the GA and the quasi-Newton method, the bulk energy given in Eq.(l 1) and Eq.(12) 
can be minimized with respect to displacements. To demonstrate the capability of this hybrid 
optimization method, three numerical experiments are reported here. For all three experiments, the 
GA first evolves a given number of generations, and the fittest chromosome in the last generation 
is selected. This fittest chromosome is then used as the initial guess for the quasi-Newton method, 
and is refined iteratively. The termination criterion for the quasi-Newton method is based on the 
assumption that the norm of difference between two consecutive potential solutions is smaller than 
the predefined value 6 (chosen in experiments as 5 = 1 x 10~^): 



(29) 
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All simulations reported here have been carried out for Au23Cu3oZn47. For this specific material, 
its physical parameters are available for the ID case: [13,25]: 

a2 = iSO g/ms^cmK, 04 = 6 x 10^ g/ms^cmK, qq = 4.5 x 10^ g/ms^cmK, 

00 = 208K, p = ll.lg/cm^, = S.mig/ms'^cmK, k = 1.9 x 10-^cmg/ms^K. 

Since 2D experimental values are not available to us, we take all the parameters in the Landau 
free energy function the same as above and complete parameterization of the model by assuming 
ai = 2a2 and 03 = 02 as suggested in [17,12]. Earlier, we confirmed numerically that the es- 
sential features of the 2D problem can be captured with this parameterization, at least in the case 
of square to rectangular transformations considered here [28,29]. Recent studies presented in [3] 
provided encouraging results also for more general cubic to tetragonal transformations. However, 
experimental results on multi-dimensional SMA samples are still lacking. The next step and a nat- 
ural development of the study presented here would be accounting systematically for the dynamics 
of the thermal field. Numerical experiments pertinent to this generalization would be much more 
involved. 

The first experiment is the simulation of phase combinations in a ID SMA wire with the bulk 
energy given by Eq.(12). The physical interpretation of the problem is sketched in Fig.2. The 
length of the wire is 1cm , and the applied mechanical force is / = 500g/ {ms^cw?) which is 
evenly distributed along the whole length. The initial condition for this case, u = (rand — 0.5)/5, 
provides a random distribution of displacement in the interval [—0.1, 0.1]. 15 nodes are used for 
the spatial discretization. For the filter operation in the GA, 7 nodes are used for the sparser grid to 
remove the higher frequency components. The GA evolves 800 generations with 60 chromosomes 
in each generation, each chromosome consists of displacement values within the range from —0.1 
to 0.1 on internal nodes. 

To analyze the performance of the hybrid optimization method for phase combination simulation, 
the bulk energy given by Eq.(12) has been monitored in the GA, and plotted in Fig.4 (left). The 
quasi-Newton iteration process has also been monitored. The actual update step size in the quasi- 
Newton method is plotted in Fig.4 (right). The trends in the two curves demonstrate how the 
optimization process evolves. As expected, the GA starts to converge to the global minimizer 
gradually, but has difficulty to locate it precisely. The quasi-Newton method starts with the output 
of the GA as its initial guess, and refined the global minimizer effectively. 

The final estimation for the phase combination in this case is given in the right column in Fig. 5: 
the order parameter e (upper plot) and the displacement (lower plot). We observe that the entire 
domain is divided into two parts, one with the strain values e ~ 0.11 and the other with e ~ —0.11. 
This phase combination agrees well with the previously reported results (e.g., [13,19,25]). For 
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comparison purpose, we have also provided an approximation obtained with the GA at the initial 

stage of the procedure (see the left column of Fig. 5). Although the quality of the final solution has 
been refined with the quasi-Newton procedure, we note that the two parts in the distribution of e 
are captured by the GA. Quantitative values are easily identified as being between 0.1 to 0.15 for 
one of the domains and between —0.15 to —0.1 for the other. 

The second experiment aims at simulating the phase combination in a 2D SMA patch sketched 
in Fig. 3. The size of the patch is [—1, 1] x [—1, l]cm^ and the applied mechanical forces fx — 
fy = SOOOg/ {ms^cm^) are distributed evenly through the entire patch. In this case, we use 
Ux — (rand — 0.5)/5 and Uy — (rand — 0.5) /5 as the initial conditions. In this experiment, 
we use different node numbers for the GA and the quasi-Newton method. In the GA, there are 9 
nodes in each direction used for the discretization along the x and y axes. This leads to the total 
number of optimization parameters being 128. The filter operation is carried out in the same way as 
explained for the ID experiment, but for both x and y directions, and the number of nodes in each 
direction on the sparse grid is 6. The GA is firstly run 1500 generations with 120 chromosomes 
in each generation. The output of the GA is then interpolated onto a denser grid with 15 nodes 
in each direction, and is used as the initial guess for the quasi-Newton method. The interpolation 
methodology is based on the pseudospectral method with the Chebyshev collocation points, as 
discussed in Section 3. 

The final phase combination estimated in this case are presented in Fig. 6. The distribution of the 
order parameter 62 is given in the upper right subplot. For comparison purpose, the estimated 
distribution of 62 from the GA is presented in the upper left subplot. The final displacement distri- 
butions along the x and y directions are presented by the two lower subplots. Once again, the entire 
domain can be divided by two parts with values e2'^Q.lle2'^ —0.11, respectively, referred to as 
martensite plus and martensite minus [2A,19\. 

In the last experiment, we consider the same patch [0, 1] x [0, l]cm^ with modified distributed 
forces fy — 0,fx — 2000g'/ {ms^cw?) (with the same initial conditions as in the previous experi- 
ment). All other computational parameters are taken the same as those in the previous experiment. 
The numerical results are presented in Fig. 7 in a way similar to already reported. In this case, the 
entire structure is divided into martensite plus and minus in a different way, the interface is the 
central vertical line due to the horizontal symmetric loading. The simulated order parameter is still 
close to either 0.11 or —0.11 and the estimated 62 distribution obtained with the GA can capture 
the essence of its final profile. 

In the analysis that follows we explain why the order parameter takes values close to either 0.11 
or —0.11. From the Landau free energy function, it is easy to estimate the order parameter value 
which minimizes the local free energy by setting: 
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For the ID case, the order parameter should be replaced by e. A simple operation gives the follow- 
ing equation: 

02 (A6')e2 - 0462 + 0562 = 0, (31) 



which means: 



62 = ^-A, ■ (32) 

lap, 



Note that 62 = is also a solution to (31) and can be associated with the austenite phase which 
is unstable at the current temperature. Therefore, martensitic phases in this case are of greater 
interest. The temperature difference from the transformation temperature here is given as = 
9 — 9q = 2", and the local minima can be estimated as: 62 = ±0.1146. 

All three numerical experiments have demonstrated that the distributions of the order parameter 
are represented by a combination of two different values, 62 = ±0.11, which agrees well with the 
prediction of the above analysis. 



5 Conclusions 



In the present paper, phase combinations in the ID and 2D SMA structures have been analyzed 
in the case of square to rectangular transformations. The phase combinations have been obtained 
by minimizing the bulk energy in the considered SMAs structures, subject to given temperature 
distribution, mechanical loadings, and boundary conditions. By combining the global and local 
search techniques, the developed hybrid optimization method provides a promising strategy for 
the solution of the associated non-convex problem. 
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Fig. 1. (a) Sketch of the square to rectangular transformation, (b) The temperature dependency of the free 
energy function for the transformation. 
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Fig. 2. Sketch of an one dimensional SMA wire under mechanical loading 




Fig. 3. Sketch of a two dimensional SMA patch under mechanical loading 
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Fig. 4. Convergence of the hybrid optimization method. Left: reduction of the bulk energy in the GA. Right: 
reduction of the step size in the quasi-Newton method 
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Fig. 6. Phase combination in a SMA patch induced by mechanical loading along the x and y directions. 
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Fig. 7. Phase combination in a SMA patch induced by mechanical loading along only the x direction 



